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ABSTRACT 

We present a simple method for fitting parametrized mass models of the Milky Way to 
observational constraints. We take a Bayesian approach which allows us to take into 
account input from photometric and kinematic data, and expectations from theoretical 
modelling. This provides us with a best-fitting model, which is a suitable starting 
point for dynamical modelling. We also determine a probability density function on 
the properties of the model, which demonstrates that the mass distribution of the 
Galaxy remains very uncertain. For our choices of parametrization and constraints, 
we find disc scale lengths of 3.00 ± 0.22 kpc and 3.29 ± 0.56 kpc for the thin and 
thick discs respectively; a Solar radius of 8.29 ± 0.16 kpc and a circular speed at the 
Sun of 239 ± 5 kms~^; a total stellar mass of 6.43 ± 0.63 x 10^" M©; a virial mass of 
1.26 ± 0.24 X 10^2 and a local dark matter density of 0.40 ± 0.04GeVcm"^ We 
find some correlations between the best-fitting parameters of our models (for example, 
between the disk scale lengths and the Solar radius), which we discuss. The chosen 
disc scale-heights are shown to have little effect on the key properties of the model. 
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1 INTRODUCTION 

■ A great deal is still unknown about the distribution of mass 
, in the various components of the Milky Way. The major dis- 

■ coveries in Galactic astronomy over the past decade have al- 
most all been related to components which comprise a small 
fraction of the total mass of the Milky Way, most of which 
either are or were dwarf galaxies (for ex ample, the many 

' objec ts observed in the "Field of Streams" , iBelokurov et al.l 
l2006h . The structure of the dominant components - the 
disc(s) and cold dark matter (CDM) halo - remains rather 
uncertain. 

An important element of understanding and constrain- 
ing the structure of the major components of the Galaxy is 
creating Galaxy models which can be compared to observa- 
tional data. It is important to draw a distinction between 
three types of Galaxy models: mass, kinematic and dynami- 
cal models. Mass models are the simplest of these, and only 
attempt to describe the density distribution of the various 
Galaxy components, and thus the Galactic potential (e.g . 
iKlvpin. Zhao, fc Somer^dU^ |2002| : iDehnen fc BinnevI 1 199^ . 
henceforth DB98). Kinematic models, s uch as those pro- 
duced by GALAXIA ISharma et al.|[201ll ). specify the den- 
sity and velocity distributions of the luminous components 
of the Galaxy, but do not consider the question of whether 
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these are consistent with a stea dy state in any Galactic po- 
tenti al. Dynamical models (e.g. IWidrow. Pvm. fc Dubinskil 
|2008| ) describe systems which are in a steady state in a given 
potential, because their distribution functions depend only 
on the integrals o f motion. The Besangon Galaxy model 
dRobin et all 1200^ ) is primarily a kinematic model with a 
dynamical element used to determine the vertical structure 
of the disc. 

It is clear that moving beyond simple kinematic models 
to full dynamical ones is an essential step in fully under- 
standing our Galaxy. The majority of the mass of the Galaxy 
is expected to lie in the CDM halo, which is only observable 
through its gravitational effect on luminous components of 
the Galaxy, so purely kinematic models cannot provide any 
insight into its structure. 

The first step towards a dynamical model is to produce 
a mass model that is consistent with available const raints. 
An influential mass model was that o f lSchmidtl(|l956l ). and, 
as observational data and understanding of galaxy structure 
improved, upd ated versions have been prod uced by other au- 
thors, notably I Caldwell fc Ostrikeij and DB98. Our 
intention in this study is to follow these authors in produc- 
ing a mass model that is consistent with up-to-date obser- 
vational data and theoretical understanding, and to provide 
a simple framework for producing these models into which 
future data can be placed as they become available. 

The major difficulty in producing a model of this kind 
is drawing together data from numerous different studies of 
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the various components which make up the Milky Way in 
a way that is as consistent as possible. Such studies often 
make different underlying assumptions, and sometimes come 
to seemingly mutually contradictory conclusions. In princi- 
ple the correct approach is to return to the raw data that 
each study was based upon and to synthesise them into one 
coherent picture. Even if this is possible in practice, it is cer- 
tainly an immense undertaking, and one we do not attempt 
here. Instead we follow the approach of previous authors by 
accepting various constraints on the parameters of our mod- 
els as stated by other studies, without returning to the raw 
data. In addition we use well understood kinematic data 
sets, and - for the CDM halo, about which observational 
data is limited - we use our best current theoretical under- 
standing. Our approach is best thought of as Bayesian with 
direct constraints on the model parameters (from photomet- 
ric data or theoretical insight) being our Bayesian "prior", 
and kinematic data used to find the likelihood. 

In this paper we present a simple method for deter- 
mining both a best-fitting parametrized mass model of the 
Galaxy, and the full probability density function of the pa- 
rameters of the model. This paper is associated with a se- 
ries of papers in wh ich we construct dynamical models to fit 
observational data (|Binnevll2O10l : iBinnev fc McMillanllioTll . 
McMillan et al. in preparation). Eventually these dynami- 
cal models will themselves be used to constrain the Galactic 
potential. 

Like DB98 we restrict ourselves to axisymmetric mod- 
els. It is clear that the Galaxy is not actually axisymmetric, 
especially the inner Ga laxy (see Sectionj2.1ll , but the disc is 
close to axisymmetric (jjuric et a l. 2008), and (as noted by 
DB98) axisymmetric models successfully account for obser- 
vations in the 21-cm line of hydrogen at Galactic longitudes 
I > 30, corresponding to J? > 4 kpc. 

To find the gravitational potential associated with a 
given mass model we use the publicly available code GALPOT, 
which is described by DB98 section 2.3. 

In Section [2] we describe the different components that 
make up the bulk of the mass of the Galaxy, and how our 
model represents them. We also explain our Bayesian priors 
on the parameters that describe our model. In Section [3] 
we give the kinematic data we use to constrain our model. 
All the constraints applied to the model are summarised in 
Table [T] In Section |4l we describe the method used to fit the 
model, and in Section [5] we give our best-fitting model and 
the rest of our results. In Section [6] we compare our results 
to other data sets. 

Unless otherwise stated, any measurement or derived 
quantity we use to constrain our model is assumed to have 
Gaussian uncertainties. 



2 COMPONENTS OF THE MILKY WAY 
2.1 The bulge 

The Galactic bulge has been shown in many studies to be 
a near-prolate, triaxial rotati ng bar with its long axis in 
the p lane of the Galaxy (e.g. iBinnev. Gerhard, fc Spergell 
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[1993). However, our interest in this study is in producing 
an axisymmetric model, so we are forced to make crude ap- 
proximations in our modelling of this component. We must 



therefore accept that it is unwise to compare our model to 
measurements taken from the inner few kpc of the Galaxy, 
as it can be expected to do a poor job of reproducing them. 

Insight into the structure of the bulge can be gained 
from photometric studies, however one must be careful when 
doing so as there can be a major contribution to the stel- 
lar density in the inner few kpc from the disc component. 
The model used for the disc in these studies will there- 
fore have a significant effect on the properties determined 
for the bulge. This goes some way toward s explaining why 
the mass of the bulge as determined by IPicaud fc RobinI 
(2004), 2.4 ±0.6 X 10^° M©, using a photometric model with 
a "hole" in the centre of the disc, is so much larger than 
that determined by studies using kinematic data (e.g. pB98 ; 
IWidrow. Pvm. fc Dubinskill2008l : Isissantz fc Gerhardll2002h 
which do not include a disc "hole" in their models. 

O ur density profile is based on the parametric model 
which iBissantz fc Gerhard! (|2003) fit to dereddened L- band 



COBE/DIRBE data (ISpergel. Malhotra. fc Blitl 19961) . and 
the mass-to-light ratio that IBissantz fc GerhardI determine 



from a comparison between gas dynamics in models and 
those observed in the inner Galaxy. This model is found on 
the assumption that the disc component has no central hole. 
This also assumes that the mass-to-light ratio is spatially 
constant, which allows us to convert a photometric model 
directly into a mass model for this component. 

The IBissant z fc Gerhard model is not axisymmetric, so 
we make an axisymmetrised approximation which has the 
density profile 



Ph = 



Ph,0 



exp 



(l + r'/ro)'' 
where, in cylindrical coordinates, 



(1) 



(2) 



with a = 1.8 , rp = 0.075kpc, rcut = 2.1kpc, and axis ratio 
q = 0.5. The IBissantz fc GerhardI mass-to-light ratio has a 
quoted uncertainty ±5%, but given the extent to which we 
have altered this model it is prudent to recognize the need to 
introduce further uncertainty in our model fitting. We there- 
fore assume that the bulge mass Mb = 8.9 x 10^ Mq, with 
uncertainty ±10%. For this density profile, this corresponds 
to scale density /9b,o = 9.93 x 10^" Mq kpc"^ ± 10%. 

2.2 The disc 

The Milky Way's disc is usually considered to have two 
major components: a thin disc and a thick disc (e.g. 
iGilmore fc Reidlll983D . These are generally modelled as ex- 
ponential in the sense that 



Pd(R,z) = - — exp — 

2zd V -Zd Hd 



(3) 



with scale height z^, scale length Rd and central surface 
density Ed,o. The total mass of a disc like this is Md = 
27rEd,o^?d- 

As mentioned in Section 12.11 some studies of the in- 
ner Galaxy prefer models w ith a central "hole" in the disc 
(e.g. iPicaud fc Robinll2004h . We do not consider such mod- 
els here. The kinematic data we consider (Section |3]) is all 
related to parts of the Galaxy which lie outside any cen- 
tral disc hole, so, in this study, a model with a central hole 
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should simply redistribute mass in the inner few kpc from 
the disc to the bulge (note that our prior on the bulge would 
have to be replaced, as it is taken from a study which used 
a model of the disc with n o central hole). 

The lJuric et alj (|2008l ) analysis of data from t he Sloan 
Digital Sky Survey fSDSS: lAbazaiian et alll2009l ) showed 
that the approximation to exponential profiles is a sensible 
one for the Milky Way, and produced estimates based on 
photometry for the scale lengths, scale heights and relative 
densities of the two discs. 

The scale-heights of the discs are not at all well con- 
strained by the kinematic data we use in this study, so 
initially we accept without question the lJuric et al.l (|2008l ) 
best-fitting values, Zd.thin ~ 300 pc and 2d,thick = 900 pc. In 
Section FS.ll we explore the (relatively small) effects of chang- 
ing the assu med disc scale-heights. 

The iJuric et al.l scale-lengths for the thin and thick 
discs are 2.6 and 3.6 kpc respectively, with a quoted un- 
certainty of 20% in each case. The local density normalisa- 
tion /d,0 = pthick(-R0,20)/pthin(-R0,20) is quoted as 0.12, 
with uncertainty 10%. We approximate these uncertainties 
as Gaussian and uncorrelated, and take these as prior prob- 
ability distributions on i?d,thin, -Rd.thick and /d,0- Again, we 
assume a constant mass-to-light ratio which allows us to con- 
vert these photometric constraints directly into constraints 
on the mass density. 



2.3 The dark-matter halo 

For obvious reasons, we cannot use photometric data to con- 
strain the shape of the dark matter distribution, so instead 
we look to cosmological simulations for insight. 

In cosmological simulations that only include dark mat- 
ter, halo density profiles ar e well fit by a universal profile , 
known as the NFW profile (|Navarro. Frenk. fc Whitdll996l ') 



Ph = 



Ph.O 



•.{i + xy 



(4) 



where x = r/rh, with Vh the scale radius. It is clear that 
the baryons within CDM haloes will have an impact on the 
halo profile, but the nature of this impact remains uncer- 
tain because it depends on the complex physics of baryons. 
Cosmological simulations suggest that the condensation of 
baryons to the centre of galactic haloes will cause the den- 
sity profile in the inner parts of t he halo to be stee per than 
p oc r"\ as it is in the NFW case jPuffv et al.ll2010l i. Obser- 
vations of low surface brightness galaxies, however, appear 
to indicate that they lie in dark-matter haloes with constant 
density cores. This difference between observation and sim- 
ulations (with or without baryonic physic s) is known as the 
"cusp-core problem" (for an overview see lde Blokll201o l). 

Observations of low surface brightness galaxies and 
dwarf galaxies can used to provide more reliable information 
about the inner dark matter profile than galaxies like the 
Milky Way because they are dark matter dominated right 
to their centre, so the mass contribution of the baryons can 
be modelled with little uncertainty compared to the total 
mass. The Milky Way is not dark matter dominated and it 
is very hard to determine slope of the dark matter density 
profile in the inner Galaxy because the baryonic density is 
dominant in the inner Galaxy, and still uncertain. 



Haloes in dark-matter-only cosmological simulations 
tend to be significantly pr olate, but with a g reat deal of 
variation in axis ratios fe.g. lAllgood et ai]|2006l ). It is recog- 
nised that, again, baryonic physics will play an important 
role - condensation of baryons to the centres of haloes is 
expected to make them rounder than dark-mat ter-only sim- 
ulations would suggest iDebattista et af]|2008l ). The shape 
of the Milky Way's halo is still very much the subject of 
debate, with different efforts to fit models of the Sagittar- 
ius dwarf's orbit favouring conflic ting halo shapes (see e.g. 
iLaw. Maiewski. fc Johnstonllioogj ). 

In this study we do not attempt to address the ques- 
tion of the shape of the Milky Way halo, or the impact of 
baryonic physics on the CDM density profile. We take the 
simple model of a spherically symmetric NFW halo (eq. |4]). 
As better understanding of the effect of baryonic physics on 
CDM haloes or constraints from observations become avail- 
able, these simple assumptions will have to be revisited. 

It is common to describe CDM haloes in terms of their 
virial mass Mvir - which is the mass contained within the 
virial radius rvir - and the concentration c — r_2/rvir, 
where r_2 is the radius at which the logarithmic slope of 
the density profile dlogp/dlogr = —2 (for an NFW pro- 
file, r_2 = rh). The virial radius Tvir is defined as the radius 
of a sphere centred on the halo centre which has an average 
density of A times the critical density, but the definition of 
A varies between authors. The definitions that are relevant 
to this paper are A — 200 (when using this definition we will 
use the notation r„ and Cv for the virial mass, radius, 
and the concentration), and a value of A which varies with 
redshift, derived from spherical top-hat collapse models (e.g. 
iGunn fc Gott .1972 ). with A ~ 94 today (for which we use 
the notation M „ /, r^i , c„/). 

IGuo et al] (|2010l ) compared the distribution of halo 
virial masses f ound in the Millenn i um and Millennium- 
II simulations (|Springel et al.l l2005l : IBovI an-Kolchin et al.l 
|2009|) to the d is tribut ion of galaxy stellar masses found 



bv iLi fc Whit^ (l2009h usi n ;^ Sloan D i gital Sky Survey 
aJl2 



data jAbazaiian et al.l 120091 ). ICuo et all (|2010l ') found that 
the ratio of total stellar mass M* to halo mass M„ was well 
fit by 



Mo 



+ 



Mo 



(5) 



with an intrinsic scatter of order 0.2 in logj^Q M«, where A = 
0.129, Mo = lO^^ "* M0, a = 0.926, /3 = 0.261 and 7 = 2.440. 

To gain insight into the likely val ue of the concentration 
parameter c„/ we consider the study of lBovlan-Kolchin et al.l 
(2010). This examined haloes taken from the Millennium- 
II simulations at redshift zero, in the mass range 

- a mass range that the Milky Way's 
halo is likely to lie in - and determined that the probability 
distribution of the concentration was well fit by a Gaussian 
distribution in In c„' , with 



Incj,/ = 2.256 ±0.272. 



(6) 



Cosmological simulations predict that typical halo concen- 
tration varie s with halo mass, but typically only weakly (e.g. 
Ci, oc M.r^ '^. lNeto et al]|2007l ) when compared to the intrin- 
sic scatter in c at a given virial mass, so we acc ept the re- 
lationship given bv iBovlan-Kolchin et all l|2010l ) as stated. 
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independent of M„/ . Baryonic physics is likely to have an 
effect on the c oncentration, much as it does on the inner 
density profile ijPuffv et al.ll20f(t l. but we neelect that in 
this study. 

The Galaxy's dark-matter halo is far more m assive than 
its st ellar halo, which has a mass ~ 4 x 10* Mq (jBell et alJ 
l2008h . We therefore treat the stellar halo as a negligible frac- 
tion of the Galactic halo, and do not consider it further. 



3 KINEMATIC DATA 

3.1 The Solar position and velocity 

If we do not know the position and velocity of the Sun it 
is extremely difficult to interpret any other kinematic ob- 
servations of the Milky Way. Unfortunately there is still 
significant uncertainty on the key parameters, namely the 
distance from the Sun to the Galactic Centre _Ro, the ro- 
tational speed of the local standard of rest (LSR) no, and 
the velocity of the Sun wi th respect to the LSR v© (see e.g. 
iMcMillan fc Binnevllioiol ) . In this work we take Ro from the 
study of stellar orbits around the supermassive black hole 
at the Galactic Centre, Sgr A*, bv iGillessen et al.. (|2009|) 



Ro = 8.33 ± 0.35 kpc. 



(7) 



thou gh we could equally have taken the value 8.4 ± 0.4 kpc 
from lGhez et al.l (|2008h . found from similar data. 

(iBinnev 



A series of recent 



pape 



2010l: 



McMillan fc BinnevI l2010l : ISchonrich. Binnev. fc DehnenI 
2010h have re-examined evidence regarding the 
value of vq. We use the value determined by 
ISchonrich. Binnev. fc DehnenI : 



V0 = {Uq,Vq,Wq) 

= (11.1, 12.24, 7.25) ± (1, 2, 0.5) kms"^ 



(8) 



where Uq is the velocity towards the Galactic Centre, Vq is 
the velocity in the direction of Galactic rotation, Wq is the 
velocity perpendicular to the Galactic plane. We quote the 
systematic uncertainties here, because these dominate the 
total uncertainty (especially in Vq). 

Many studies that can be used to constrain the value 
of vo are ones which actually constrain the ratio vo/Ro ~ 
for example using the Oort constants A and B, which can 
be determined locally (e.g. iFeast fc Whitelo ck 1997) where 
A ~ B — vo/Ro- In this study, we use the proper motion 
of Sgr A* in the plane of the Galaxy, as determined by 
iReid fc Brunthale j (|2004l l 



MSgrA* = —6.379 ± 0.026 mas yr ^. 



(9) 



Since Sgr A* is expected to be fixed at the Galactic Centre 
to within ~ lkms~^, this proper motion is thought to be 
almost entirely due to the motion of the Sun around the 
Galactic Centre {vq + Vq)/Ro. This measurement is suffi- 
ciently accurate (and the assumed velocity of Sgr A* suf- 
ficiently small) that the greatest uncertainty in the value 
of Vo/Ro required for our models actually comes from the 
uncertainly in the value of Vq ! 

This measured prope r motion is inconsistent w ith the 
value of A - B found by iFeast fc Whitelockl ([1993) - who 
give the most commonly used values for A and B. Therefore 
we do not attempt to use the Oort constants as constraints. 



3.2 Terminal velocity curves 

For circular orbits in an axisymmetric potential, the peak 
velocity of the interstellar medium (ISM) along any line-of- 
sight at Galactic coordinates 6 = for —90 < I < 90 cor- 
responds to the gas at the tangent point, at Galactocentric 
radius R — Rosinl. This is known as the terminal velocity, 
and is given by 



Vte 



VciRo sinZ) — i'c(-Ro) sinZ. 



(10) 



iMalhotral l|l994l . Il995l ;) produced data which gave this peak 
velocity for a range of lines of sight in the Galaxy from 
a number o f surveys jWeaver fc Williamd 1974 Kerr et al] 



ll986l : lBania fc Lockmanll 19841 : iKnapp et al.lll9'85F K 

To constrain our model, we have to take account of 
the effects that non-axisymmetric structure in the Galaxy 
and non-circular motion of the ISM will have on these data. 
To do so we follow DB98 in allowing for an uncertainty of 
7kms~^ in «tcrm at any given I, and restricting ourselves 
to data at | sinZ| > 0.5, which is unlikely to be significantly 
affected by the bar. 



3.3 Maser observations 

In recent years it has become possible to use Galactic 
maser sources as targets for astrometric measurements that 
are sufficiently accurate for parallaxes with uncertaintie s 
~ lO/xas to be determined for them (e.g. iReid et "al]|2009l ). 
iMcMillan fc BinnevI l|2010l) showed that these observations 
were consistent with models that placed the masers on near 
circular orbits with vo / Ro similar to the value implied by the 
proper motion of Sgr A* (Section |37l}. However Ro by itself 
(or, equivalently, vo) was shown to depend quite strongly 
on the shape of the rotation curve. We constrain our mass 
model using these dat a by applying a vers i on of the likeli- 
hood analysis used bv IMcMillan fc BinnevI (|20ld ). 

The likelihood analysis for these data requires integrat- 
ing over heliocentric radius the probability density of the 
observations (parallax, proper motion and radial velocity) 
given a model for the expected velocity distribution. The 
maser sources are associated with high mass star forma- 
tion regions, and are expected to be on near-circular or- 
bits. We model the velocity distribution as circular rota- 
tion (with a velocity which depends on the potential), plus 
a random component that has a Gaussian probability dis- 
tribution, uniform in direction, with width A^,. The inter - 
ested reader can find details in IMcMillan fc BinnevI (|20lO). 
However the analys i s in this study differs from that of 
IMcMillan fc BinnevI l|2010l) in that (a) the rotation curve, 
including no, is defined by the mass model; (b) the value 
of V is assumed to be that in eq. ( c) we include data 
from iRvgl et all 1I2OIOI) and ISato et al.1 (I2OI0I). w hich were 
not available when the lMcMillan fc Binnev ( 2010l ) study was 
performed; (d) in the interests of simplicity we fix A„ = 
7kms~'^, close to the best-fitting values found when A„ was 
allowed to vary in the previous study, and (e) we assume the 
masers have no systematic velo city offset from the rotation 
curve - one of the conclusions of IMcMillan fc BinnevI (|2010D 
was that no such offset is required to explain the data. 
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3.4 Vertical force 

iKuiiken fc Gilmorel (|l99ll ) used observations of K stars to 
find tlie vertical force at 1.1 kpc above the plane at the Solar 
radius, -fT^^i.i. We adopt the value they found 



ifz.i.i = 27rG X (71 ± 6) Mo pe- 
as a constraint. 



(11) 



3.5 Observational constraints on the structure at 
large radii 

It is extremely difficult to gain useful constraints on the 
structure of the Milky Way at large radii from observational 
data. Any population of dynamical tracers suffers from small 
number statistics and/or poor observational constraints (es- 
pecially on proper motion) and because of uncertainty over 
the associated distribution function, especially the veloc- 
ity anisotropy, and whether or not all of the population is 
in fact bound to the Milky Way. This leads to total halo 
mass estimates which have fractional uncertainties of o rder 
unity, e.g. 1.9t? j x 10^^ Mo (Wilkinson fc Evansll 19991 ). 0.3 
to 2.5 X 10^^ Mq (|Battagha et al.ll2006l '). 

It is possible to take cosmologically motivated simula- 
tions and use them to provide some insight into th e expected 
distrib ution function of a tracer population. IXue et ali 
l|2008l) did this for a sample of blue horizontal branch (BHB) 
stars and found that the mass at radii less than 60 kpc, under 
this assumption, was M(< 60 kpc) = (4.0 ± 0.7) x lO" Mq. 
However this value is entirely dependent on the assumption 
that these cosmological simulations (and their prescriptions 
for star-formation, feedback etc.) can be relied upon to pre- 
dict accurately the distribution function of BHB stars. 

Another possible source of information is the motion of 
the Magellanic clouds, but this can only be used straight- 
forwar dly if they c a n be assumed to be bound to the Milky 
Way - iBesla et all (|2007l ) have shown that this may not be 
a valid assumption. 

Other authors (notably ISmith et al.l |2007| ) have at- 
tempted to determine the local escape velocity from the 
Milky Way, and use this to constrain the mass distribution 
at large radii. However this is determined under the assump- 
tion that the population of stars close to the escape velocity 
can be described by a steady-state distribution function that 
is depends on energy alone. Stars that have velocities that lie 
close to the escape velocity are on orbits with very long pe- 
riods, which means that there is no reason to expect them 
to be in a steady state, especially since it is likely that a 
large fraction of halo stars are associated with streams. The 
assumed distribution function is isotropic (because it only 
depends on energy), and this assumption is unlikely to be 
valid - halo stars in cosmological simulations typically have 
a distribution function which is radially anisotropic. 

The "Timing Argument" has been used to constrain the 
mass of the local group or that of the Milky Way itself. In its 
simplest form, this assumes that the currently observed ra- 
dial velocity of M31 towards the Milky Way must be due to 
the gravitational interaction of the two galaxies overcoming 
the overall cosmic expansion, and that the corresponding 
orbital motion must have happened within the age of the 
universe, and shows that this implies a low er limit on the 
mass of the two galaxies. iLi fc Whitel l|2008l ) extended this 



argument by searching cosmological simulations for compa- 
rable galaxy pairs that were moving towards one another in 
the simulations. They also looked for galaxy pairs compa- 
rable to the Milky Way and Leo I, which are moving away 
from one another, and used the Timing Argument under 
the assumption that Leo I is moving towards apocentre for 
a second time. This yielded a probability distribution on M„ 
with lower and upper quartiles [1.78, 3.09] x 10^^ Mq and a 
95% confidence lower limit of 7.94 x 10^^ Mq. This approach 
rather relies upon Leo I being bound to the Milky Way, 
which it may not be, and the authors note that, compared 
to the Milky Way/M31 system, the "complex dynamical sit- 
uation offe rs greater scope for u ncert ainty". 

While IWilkinson fc Evand (|l999l ) quoted a very uncer- 
tain figure for the total halo mass, they also found that the 
mass within 50 kpc was a more robustly determined quan- 
tity: A'/so = 5.4l3;2 X 1O"M0. The quoted uncertainty is 
very asymmetrically distributed about the most likely value, 
and we choose to see this result as providing an upper bound 
on the mass inside 50 kpc, which we take into account in our 
analysis through the probability distribution 

for M50 ^ MwE 

for A/50 > MwE ^^"^^ 



where AfwE = 5.4 x lO" Mq, Sm^o = 2 x Mq, and C 
is a normalisation constant. The only effect of using this 
P{Mz,o) constraint is to penalise models with Afso > A/we. 
For models with Afso <C AfwE this would be an inappropri- 
ate choice but, as we show in Section (5] that is not the case 
for this study. 

Given the difficulties described here, we have decided 
not to constrain our model with the results of any other 
study. We do compare our best-fitting model to them (Sec- 
tion [6} 



4 FITTING THE MODELS 

We use Bayesian statistics to find the probability density 
function (pdf ) for our model parameters given the kinematic 
data described in Section O and the prior probabilities given 
in Section [2] (and the value of _Ro, given in Section [3)|. We 
refer to the parameters collectively as 0, and the data as d. 
Bayes theorem tells us that this pdf is then 



p{0\d) = 



c{d\e)p{0) 
p{d) 



(13) 



where the total likelihood jC{d\0) is the product of the likeli- 
hoods associated with each kinematic data-set or constraint 
described in Section [31 and p{0) is the probability of a 
parameter set given the prior probability distributions de- 
scribed in Section [2] (and that on Rq). The Bayesian Evi- 
dence, p{d), is often a very important quantity, but in this 
study it is an unimportant normalisation constant, so we 
ignore it. 

Given the constraints we have applied to the compo- 
nents described in Section [2l we are left with 8 model 
parameters that we allow to vary: The scale-lengths 
and density normalisations of the thin and thick discs 
(-Rd,thin, Sd.o.thin, /?d,thick, Sd,o,thick); the density normalisa- 
tion - and thus mass - of the bulge (pb,o); the scale-length 
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Property constrained Constraint 


Section described 


Source 


Bulge profile 

Mb 
Disc profile 

^d.thin 
^d, thick 
-^d,thin 
^d, thick 
/d,0 
Halo profile 
M./M„ 
In c„/ 
Ro 

A^SgrA* 


see equation |[TJ 
(8.9 ±0.89) X IO^Mq 
Double exponential 
0.3 kpc 
0.9 kpc 

3.6 ±0.72 kpc 
0.12 ±0.012 
NFW profile 
see equation JSj 
2.256 ±0.272 
8.33 ± 0.35 kpc 
-6.379 ± 0.026 mas yr-l 
2ttG X (71 ± 6) Mq pc-2 
< 5.4 X 10" Mq, see eq. |[T2t 


[2:21 
[2:21 
[22] 
[2:21 
[221 
[221 
[Ql 

[231 

[273I 

ixn 

[331 
1331 


Bissantz & Gerhard ('2002') 
Bissantz & Gerhard ('20021 

Juric et al. ('2008) 
Juric et al. (2008") 
Juric et al. (2008) 
Juric et al. (20081 
Juric et al. (2008) 
^avarro et al. ( 1996) 
Li & White (2009) 
Bovlan-Kolchin et al. (20101 
Gillessen et al. (20091 
Reid & Brunthaler (20041 
Kuiiken & Gilmore (19911 
Wilkinson & Evans (1999) 




Kinematic data 


Section described 


Source 




Terminal velocities 
Maser observations 


[321 
[331 


Malhotra ('1994. 19951 
Reid et al. (2009): Rvel et al. (2010): 
Sato et al. (2010) 



Table 1. Summary of all the constraints applied to the models. 



and density normalisation of the CDM halo (rh,ph,o); and 
the solar radius (Ro). While each of these parameters is free 
to vary, each one is explicitly or implicitly associated with 
at least one prior probability distribution. 

To explore the pdf p(0|d) we use the Metropolis algo- 
rithm (|Metropolis et ah! Il953l l , which is a Markov Chain 
Monte Carlo method for drawing a representative sample 
from a probability distribution. This allows us both to find 
the peak of the pdf (to reasonable accuracy) and to charac- 
terise its shape. We start with some choice for the parame- 
ters On, and calculate p(0„|d). We then 

(i) choose a trial parameter set G' by moving from 0„ in 
all directions in parameter space, by an amount chosen at 
random from a "proposal density" Q{0' ,9n) (see below); 

(ii) determine p(0' I d); 

(iii) choose a random variable r from a uniform distribu- 
tion in the range [0,1]; 

(iv) ii p{0'\d) /p{6n\d) > r, accept the trial parameter set, 
and set On+i = 0' . Otherwise set On+\ = On- 

(v) Return to step (i), replacing 0„ with On+i- 

The first few values of are ignored as "burn-in", which 
helps to remove the dependence on the initial value of 0. 

The pdf is quite narrow in some directions in the pa- 
rameter space, but these directions are not simply paral- 
lel to the coordinate axes. Therefore it is efficient to use a 
proposal density which is aligned (to a reasonable approx- 
imation) with the pdf - this is allowed by the Metropolis 
algorithm as the only constraints on Q(0i,02) are that it 
is symmetrical with respect to swapping 0\ and O2, and 
makes it possible to reach any point in phase space. We 
therefore use the Metropolis algorithm in two phases - first 
with a proposal density which is simply made up of Gaus- 
sian distributions in each parameter individually, with asso- 
ciated step size chosen by hand. The resultant chain (of a 
relatively short length, and after a burn-in) is then used to 
construct a covariance matrix that is then used to define a 
new Q{0' ,0n). This new proposal density is, again, a multi- 
variate Gaussian but with the principal axes now aligned 



to the eigenvectors of the covariance matrix, and with the 
step size in each direction related to the respective eigenval- 
ues. The chain produced with this second proposal density 
is then used in all calculations. 



5 RESULTS 

In Figures[T]to|3]we plot the probability density functions of 
various quantities associated with our models, marginalised 
over all parameters, as determined by the Metropolis algo- 
rithm. Where appropriate we also plot the prior probability 
distribution directly associated with each. The value asso- 
ciated with our best-fitting model is indicated on each plot 
with a dashed vertical line. In Table [2] we give the parame- 
ters of our best-fitting model, and those of what we call the 
"convenient" model (see below), along with the mean and 
standard deviation for each parameter, marginalised over 
other parameters in the pdf. 

The marginalised distributions do not give a sense of 
the correlations between parameters, so in Table [3] we show 
the correlation matrix of the various parameters 0. For the 
i,jth component, corr{9i,9j), this takes the value 

, , cov(6'i, Oj) 
coTT(6i,9j) = ^ — (14) 

where cov{9i, 9j) is the covariance. A value of 1 corresponds 
to a perfect correlation, —1 corresponds to a perfect anti- 
correlation, and corresponds to no correlation. The corre- 
lation matrix is manifestly symmetric, so in Table [3] we only 
show half of it. 

The strongest correlations or anti-correlations are be- 
tween parameter pairs that describe a single component - 
this explains, for example, why the spread in M* is so much 
smaller than the spread in Ed, 0, thin or Ed, 0, thick, and why 
the standard deviation of ph,o can be as large as ~ 50% 
while that in Mv is 20% (and that in M50 is even smaller, 
see below). Other fairly strong correlations are noticeable 
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^d.O.thin 


-Rd.thin 


Sd,0, thick 


-Rd, thick 


Pb,0 


Ph.O 


''h 


Ho 


Best 




OiD.D 


o on 


zuy.o 


o.oi 


yo.D 


o noQ/ifi 
U.UUo4D 


on o 


c oo 


Convenient 




753.0 


3 


182.0 


3.5 


94.1 


0.0125 


17 


8.5 


Mean 




741 


3.00 


238 


3.29 


95.5 


0.012 


18.0 


8.29 


Std. Dev. 




123 


0.22 


110 


0.56 


6.9 


0.006 


4.3 


0.16 




DO 


Mb 






A'^.i.i 


Sd,0 


P*,© 


Ph.e 


/d,0 9t,@/9h,e 


Best 


239.1 


8.97 


66.1 


1400 


77.7 


63.9 


0.087 


0.0104 


0.122 1.38 


Convenient 


244.5 


8.84 


65.1 


1340 


75.4 


60.3 


0.083 


0.0111 


0.121 1.12 


Mean 


239.2 


8.96 


64.3 


1260 


76.5 


62.0 


0.085 


0.0106 


0.120 1.29 


Std. Dev. 


4.8 


0.65 


6.3 


240 


5.3 


7.6 


0.010 


0.0010 


0.012 0.30 



Table 2. Parameters (upper) and derived properties (lower) of our best-fitting model, our "convenient" model, and mean and standard 
deviation marginalised over all the (other) parameters in the pdf. The density profile of the bulge follows the description in eq. JlJ, 
with other parameters as described directly below that equation; the density profiles of the discs follow eq. l(3]l, with scale-heights 0.3 
and 0.9 kpc for thin and thick discs respectively; and the halo density profile is that given in eq. Sd,0 is the disc surface density 
at the Sun; p*,0 and q are the local densities of the stellar component and CDM halo, respectively; Qt^Q/gh @ is the ratio of the 
gravitational force on the Sun from the stellar component (g»,0) to that from the CDM halo (gi^ q) - this is a measure of whether the 
Galaxy is disc- or halo- dominated. The means and standard deviations of the parameters are not always particularly helpful statistics, as 
they say nothing about the correlations between parameters. It should also be noted that the true uncertainties on q thim Sd,o, thick 
and p*,0 are considerably larger than the standard deviations quoted here, because the disc scale heights have been held constant (see 
Section |5.1|I . Distances are quoted in units of kpc, velocity in kms~^, masses in 10^ M0, surface densities in M0pc^^, densities in 
M0 pc"-^, and Kz^i,i in units of (27rG) X M0 pc"^. 





Sd.O.thin 


-Rd.thin 


Sd.O, thick 


-Rd.thick 


Pb,0 


Ph,0 '"h 


Sd.O.thin 


1 












-Rd.thin 


-0.727 


1 










^d,0,thick 


-0.467 


0.483 


1 








Hd,thick 


0.468 


-0.304 


-0.845 


1 






Pb,0 


-0.197 


0.167 


-0.047 


0.037 


1 




Ph,0 


-0.632 


0.328 


-0.247 


0.151 


-0.009 


1 


rh 


0.588 


-0.299 


0.221 


-0.130 


0.003 


-0.877 1 


Ro 


-0.369 


0.602 


-0.071 


0.157 


-0.017 


0.373 -0.338 



Table 3. Correlation matrix for the model parameters. A value of 1 corresponds to perfect correlation (e.g. any parameter with itself), 
— 1 corresponds to perfect anti-correlation, to no correlation. The strongest relationships are anti-correlations between parameters 
which, taken in combination, define the mass of a given component or its density at a given point - i.e. between pb ''hi between 
Sd,o,thin & -Rd,thin and between Sj^o, thick & iid.thick- 



between Ed, 0, thin and virtually every other parameter, and 
between iZd.thin and Ro- 

Figure [1] shows the pdfs associated with the disc scale- 
lengths, i?thin and i?thick. Both have peaks around 3 kpc, 
and the best-fitting values also lie quite near 3 kpc. That the 
kinematic data (primarily the terminal velocity curve) drives 
the two scale-lengths towards a similar value is unsurprising 
- the two discs have effects on the dynamics of the Galaxy 
that scale very similarly. The posterior pdf on i?thin is much 
narrower than the prior, which is not true to the same extent 
for i^thick; this is not surprising given that the thin disc is 
more massive than the thick disc so the constraints based on 
kinematic data are more likely to dominate the photometric 
constraints in the thin disc pdf. 

DB98 considered mass models with four different scale- 
lengths between 2 kpc and 3.2 kpc, and halo density profiles 
which could take on a wide range of power-law slopes in 
their inner and outer parts. They found that the models 
with disc scale-lengths of 2.4 kpc or smaller required haloes 
which, unrealistically, had density profiles which increased 
with radius (p cx: r^) in their inner parts, whereas the mod- 
els with scale-lengths 2.8 kpc or greater did not. Therefore, 
given that we require a halo profile that increases as p oc 



in its centre, we should not be surprised that this favours 
models with disc scale-lengths larger than 2.5 kpc. 

Figure [2] shows the posterior pdfs for the position of 
the Sun and the circular speed at the Sun in our models. 
The posterior pdf on Ro is much narrower than the prior, 
with a standard deviation of 0.16 kpc, compared to an uncer- 
tainty on the prior of 0.33 kpc. The posterior pdf on vq/Ru 
is of a similar width to that which would be seen if our 
only information on its value was the proper motion of Sgr 
A*, but it is slightly displaced to higher values of vo/Ro- 
This displacem ent is most likely due t o the maser data (Sec- 
tion [33J, which iMcMillan fc BinnevI (I2OI0I ') showed are best 
fit by models in which vq /Ro is slightly larger than the value 
implied by equation The corresponding value of vo, 
239.2 ± 4.8kms~^, has a s ignificantly sma l ler un certainty 
than simpl y combining the Gillessen et al.l (|2009l ) value of 
Ro and the lReid fc Brunthalerl ( 20041 ) value of psgrA* gives, 
primarily because of the tighter constraints on Ro. 

The posterior pdf on bulge masses (Figure [3l to p left) 
is close to our prior from iBissantz fc Gerhardl l|2002l ). The 
bottom right panel of Figure [3] shows that the stellar mass 
of the models is tightly peaked at values towards the high 
end of those we would expect given the virial mass, and 
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2 2.5 3 3.5 4 4.5 

Rthick (kpc) 

Figure 1. Histogram of the pdf of the thin (top) and thick 
(bottom) disc scale-lengths in our models (solid histogram) nor- 
malised over all other parameters, compared in each case to the 
prior pdf described in Section 12.21 (dotted). In each plot the value 
corresponding to our best-fitting model is shown as a dashed ver- 
tical line. In each case the posterior pdf is peaked near to 3 kpc. 



our constraint on this relationship given in equation ©. 
This can equivalently be thought of as corresponding to low 
halo masses, given the stellar masses. Compared to the prior 
pdfs taken from theoretical considerations ('Section l2.3p . the 
haloes have low masses and high concentrations (lnc„' = 
2.73 for the best-fitting model, and with mean and standard 
deviation 2.83 and 0.18 respectively). The halo is likely to be 
sub-dominant at Ro in the sense that the gravitational force 
on the Sun is mostly from the stellar component, though this 
is not the case for all models (ratio of gravitational forces 
5*, 0/311,0 = 1-29 ± 0.30). The mass inside 50 kpc of our 
best-fitting model is 5.3 x 10^^ Mq, with mean and standard 
deviation 5.1 ±0.4 x 10^^ M^. Th i s is cl ose to the mass limit 
we take from I Wilkinson fc Evanj (|l999l ). indicating that the 
latter is an important upper limit on our models. 

The local density normalisation /d,0 is 0.122 for our 
best-fitting model. The mean and standard deviation are 
0.120 ± 0.012, which is exactly what one would expect from 
our prior alone. This is due to the fact that our kinematic 
data does not provide us with any great ability to discrimi- 
nate between mass in either of the two discs. Therefore, the 
value of /d,0 almost entirely defines (with the disc scale- 
heights) our constraints on the relative contributions of the 
thick and thin disc. 



. 1 1 1 1 1 , 1 1 






1 .... 1 - 



7.5 8 8.5 9 



Ro (kpc) 



Vo/R, 



I I I I , , I I I . I . I I , I I I I 

28 28.5 29 29.5 

v„/R„ (km/s/kpc) 

Figure 2. Histogram of the pdf of Rg (top) and vq/Rq (bottom) 
normalised over all parameters (solid line). The prior pdf on Rq 
and the pdf on vq/Rq associated with the apparent motion of 
Sgr A* (Section 13.11 1 are plotted as dotted lines. As in Figure [T] 
the value corresponding to our best-fitting model is shown as a 
dashed vertical line. 

ICatena fc Ulliol 1I2OIOI ) used a rather similar method 
to the one described in this paper with the single inten- 
tion of determining the local dark matter density, finding a 
value ~ 0.39 GeV cm~'^. We have made a number of different 
assumptions, and taken into account different constraints, 
but many of the key assumptions (axisymmetry, exponen- 
tial disc and a halo profile motivated by CDM simulations) 
are the same, and we find a similar local dark matter density 
0.40 ± 0.04 GeV cm~^. It is well worth noting that the un- 
certainty we quote here is purely the statistical uncertainty 
associated with this set of parametrized models, and the 
constraints we have chosen to apply. In particular, we have 
made the approximations that the discs are well described 
by exponential profiles at all radii, and that the halo density 
profile is well described by a spherically symmetric NFW 
profile. These assumptions and others, including ones which 
are not not directly related to the dark matter profile, will 
have a significant effect on the value found for the local dark 
matter density. We expect that the systematic uncertainties 
on this value are much larger than the statistical uncertainty 
we state here. 

It is clear from Figures[T]to[3]that there is a wide range 
of models that fit our constraints almost as well as our best- 
fitting model. It is often convenient to use models in which 
certain key parameters are chosen to take simple values (for 
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Figure 4. Comparison of the terminal velocity curve predicted 
by our best-fitting model (solid curve), the "convenient" model 
(dashed red curve) and the terminal velocity measurements from 
the Galaxy (crosses). The two model curves are nearly indistin- 
guishable. 



Figure 3. Histogram of the pdf for the bulge mass M\y (top- 
left), the total stellar mass M, (top-right), the virial mass M„ 
(bottom-left) and the log of the ratio of M, to the value of M, 
predicted by equation l(5j (bottom-right). In first and last case, 
where a straightforward prior applies to the quantity shown, it is 
plotted as a dotted curve. Again, the value corresponding to our 
best-fitting model is shown as a dashed vertical line in each plot. 



example it is common to take 7?o = 8 or 8.5 kpc), and we 
provide such a model in Table [S] To take account of the cor- 
relations between parameters we chose to build this "conve- 
nient" model one step at a time. We first hold Ro constant 
at a suitable value determined from Figure [21 and find the 
pdf associated with these models, which we use to choose 
a suitable thin disc scale length i?d,thin- We then repeat 
this process to find a suitable value for i?d, thick, and then 
again for rh. Following this procedure we take Ro = 8.5 kpc, 
J?d,thin = 3 kpc, i?d, thick = 3.5 kpc and Hi = 17 kpc. 

In Figure [4] we compare the terminal velocity curve as- 
sociated with our best-fitting model, and with the "conve- 
nient" model, to the observational data (Section 13. 2p . The 
difference between the curves of the two models is very small, 
and both provide a good fit to these data. In Figure [5] we 
plot the rotation curves of the best-fitting and convenient 
models, and the decomposition into the components due to 
the bulge, discs, and halo. The two models are more clearly 
distinguished in this plot, primarily because of the larger 
value of Ro chosen for the convenient model, which results 
in a higher value for vq, because of the strong correlation 
between the two. 



5.1 Effect of changing the disc scale- heights 

In all of the models discussed thus far (and shown in Fig- 
ures[l}{5]and Tables[2]&[3| the scale- heights of the discs have 
been held constant at 300 pc and 900 pc for the thin and 
thick discs respectively. These values were cho sen because 
they are the best-fitting scale-heights given by lJuric et al.l 




R 



Figure 5. The rotation curve for our best-fitting (black) and 
convenient (red) models - note that radius is plotted logarithmi- 
cally. The solid line, labelled F, is the full rotation curve, with the 
other curves showing, in each case, the contribution of the bulge 
(B, dotted), discs (D, short-dashed) and halo (H, long-dashed). 

( |2008t ). and were held constant because we do not expect the 
kinematic data described in Section [3] to strongly constrain 
the vertical density profile. It is, however, important that we 
check this assumption by considering different scale-heights 
for the disc. Therefore we now consider models which have 
a thin disc scale-height /id, thin of 250,300 or 350 pc and a 
thick disc scale-height /id, thick of 750, 900 or 1050 pc, in all 
nine possible combinations. 

One change to our prior must be considered because of 
the effect of changing the scale-height of the disc. The local 
density normalisation, /d,© = pthickiRe, zq)/ pthin{R@, zq), 
is defined at the Sun's position, so if the scale-heights 
change, the value of /d,0 changes without the surface den- 
sities of the two discs changing. It has been noted (by, 
for example, iRevle fc Robiiil I2OOII ) that there is an anti- 
correlation between the values of /d,0 and /id. thick found 
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by studies using star counts, which can indeed be seen in 
figures 21 & 24 of ljuric et all l|2008l ). We therefore approxi- 
mate that /d,0 simply changes with /id, thick, taking the val- 
ues /d,0 = 0.14, 0.12, 0.10 for /id.thick = 750, 900, 1050 pc re- 
spectively, with an uncertainty in /d,© of ±0.012 in all cases 
(a s previously ) . The se values are taken by eye from figure 21 
ofQ uric et al] (|2008l ). with a correction for the recognised bi- 
ases in the data. We ignore any possible correlation between 
/d,0 and /id, thin - an anti-corre l ation of this type is seen in 
figures 21 & 24 of I Juric et al.l (|2008h . but could easily be 
related to the relationship between /d,© and /id, thick, as the 
values they derive for /id, thin and /id, thick are also correlated. 

We give the mean and standard deviation for the pa- 
rameters (and derived quantities) of all these models in the 
appendix. Here we discuss the most significant findings. 

Most importantly, the changes in scale heights have very 
little impact on the overall structure of the Galaxy models. 
The bulge mass and virial mass of the Galaxy are virtually 
unchanged, and the total stellar mass changes by at most 
~ 5%, which is much less than the standard deviation of 
the value. The disc scale lengths are also barely changed 
(changes of at most ~ 2% in i?d,thin and ~ 4% in i?d, thick). 
The surface density at the Sun is also largely unchanged, 
with changes of ~ 1% in Kz,i.i and the local disc surface 
density Ed,© changing by ~ 5%. 

The only significant change is in the values of Ed, o, thin 
and Ed, 0, thick, which show a transfer of mass from the thick 
disc to the thin disc as /id.thin increases. The mean value of 
Ed, 0, thin rises from ~ 700 M© pc~'^ for models with /id, thin = 
0.25 to ~ 800 M© pc^^ for those with /id, thin ~ 0.35, while 
the value of Ed,o,thick declines from ~ 280 M© pc~^ to 
~ 210M©pc~^ over the same range. This then affects the 
fraction of the total disc mass found in each disc, and the 
value of /d,©. 

In every case the posterior pdf of /d,© is almost precisely 
the same as the prior. This indicates that the prior on /d,© 
is, in all cases, the dominant factor determining the how the 
Galaxy model's total disc mass is divided between the two 
discs. 

We can therefore see that holding the disc scale-heights 
at fixed values to produce the statistics in Table [2] does not 
significantly alter the quoted standard deviations, except for 
parameters related to the division of the disc material into 
the thick and thin discs. Here the true uncertainies are sig- 
nificantly larger than the quoted values because the uncer- 
tainty in scale-heights must be taken into account. 



6 COMPARISON TO OTHER STUDIES 

The thin disc scale le ngth of our best-fi tting model is at 
the larger end of the lJuric et al.l (|2008l ) range, which we 
used as a prior. It is therefore also larger than the values 
in t he range 2 to 2.5 kpc found by some other recent stud - 
ies (|Oiha et al.lll999l : IChen et al.lll999l : [Siegel et al.|[2002D . 
and the value used for the old (age > 0.15 Gyr) disc i n 
the Besangon Galaxy model (2.53 kpc. iRobin et al.ll2003l ). 
On the other hand it lie s close to the scale length found 
by iKent. Dame, fc Faziol |l993), 3 ± 0.5 kpc, and towards 
the lo wer end of the range found bv lLopez-Corredoira et akl 
(|2002l ) 3.3l°;4kpc. The stellar mass of our best-fitting 
model is 6.61 x 10^° M© (with mean and standard de- 



viation 6.43 ± 0.63 X 10^"). This compares to the value 
6.1 ± 0.5 X 10^° M© found a s a "back of the envelope" esti- 
mate bv lFlvnn et al.l (|2006l ). The local disc surface density 
(Ed,© = 63.9 M© pc"^ best-fitting; 62.0 ± 7.6 M© pc"^ mean 
± standard deviation) is somewhat larger than suggested by 
studies which count i dentified matter in the Solar neighbour- 
hood, such as that of lKuiiken fc Gilmorel lll989h who found 

Flvnn et"aiT (|2006l ) who found 



Ed,© = 48 ± 8M©pc 
Ed©~49M ""-2 



1© pc 



Our best-fitting bulge mas s is close to the centr e of th e 
prior distribution we take from lSissantz fc GerhardI (|2002h . 
This is significantly larger than most of the values found 
by DB98, and comp arable to the bulge masses found by 
IWidrow et al.l l|2008l ). This is somewhat smaller than the 
masses determine d by many purely p hotometric studies ( e.g. 
Dwek et al.l 1 19951 : iPicaud fc RobinI [20041 : iLaunhardt et all 
2002h . but, as noted in Section 12.11 there is an important 



dependence on assumptions made about the disc. We note 
that the stellar mass within the inner 3 kpc in our best- 
fitting model, ~ 2.4 x 10^° M©, is ver y close to the bulge 
mass found by iPicaud fc RobinI (|2004l ) assuming a central 
hole in the disc: 2.4 ± 0.6 x 10^° M©. 

The baryonic Tully-Fisher relation describes the rela- 
tionship between total baryonic mass Mbaryon and the cir- 
cular speed Vc a t some radius, observed in external galaxies. 
iMcGaugh et al.l l|2010l ) found log^p Mbaryon = 4.0 log^o K -I- 
1.65 for disc galaxies, with a scatter that is entirely 
consistent with observational uncertainty, where 1.1 x 14 
was the actual observed circular speed. Our models have 
login A^baryon " (4.0 log^,, Vc + 1.65) = -0.19 ±0.05, suggest- 
ing that the Milky Way has a significantly higher rotational 
speed (or, equivalently, lower baryonic mass) than the Tully- 
Fisher relation predicts. A similar offset from the baryonic 
Tully-Fisher relationship was found for the Milky Way by 
iFlvnn et al.l (|2006l ). 

In Section 13.51 we noted the difficulty of finding rigor- 
ous constraints on Galactic structure at large radii, and de- 
scribed some of problems with commonly cited constraints. 
It is still instructive to compare our models to the results 
of these previous studies. The escape speed at the Sun of 
our best-fitting model is 622kms~^, with mean and stan- 
dard deviation over all models of 606kms~^ and 26kms~^ 
respectively. This compares to the quote d 90% confidence 
interval 498 to 608 kms'^ of ISmith et"all (|2007.) . The mass 
inside 60 kpc is 6.2 x lO'^^ M© (best-fitting), with mean and 
standard deviation 5 . 9 ±0.5 x 10^^ M©. This is rather larger 
than the IXue et al.l (|2008l ) value of 4.0 ± 0.7 x 10" M©. 
IXue et aL used cosmological simulations to predict the ve- 
locity anisotropy of the BHB population they were studying 
- this leads to predictions of rather high radial anisotropy, 
which drives the mass estimate towards lower values than 
more isotropic or tangentially anisotropic velocity distribu- 
tions would suggest. DB98 adopted the constraint on Mioo, 
the mass inside 100 kpc, Mioo = 7 ± 2.5 x 10"M©, based 
on then-available data. Our best-fitting model (A/ioo = 
9.0 x 10^^ M©) and the ensemble of models as a whole (mean 
± standard deviation: 8.4 ± 0.9 x 10^^) fit well within this 
range. The virial mass of our best-fitting model M„ = 1.40 x 
10^^ M© (mean ± standard deviation: 1. 26±0.24 x 10^^ M©) 
lies w ell within the range s quot ed by IWilkinson fc Evanj 
(|l999l ) or iBattaglia et al.l (|2006h . and slightly below the 
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lower quartile of the pdf given bv lLi fc Whit3 diooi ). but 
well above their 95% confidence limit. 



7 CONCLUSIONS 

We have presented a simple Bayesian method for apply- 
ing photometrically and kinematically derived constraints, 
and theoretical understanding of galaxy structure, to 
parametrized axisymmetric models of the Galaxy in order 
to investigate the distribution of mass in its various com- 
ponents. We have applied this method to models with an 
axisymmetric bulge, exponential discs and an NFW halo. 

The method we have described is sufficiently general 
that it could be applied to any sensible parametrized ax- 
isymmetric mass model, with a wide range of kinematic, 
photometric, theoretical or other constraints. The specific 
constraints we apply are summarised in Table [1] We have 
shown that these constraints still allow a wide range of 
Galaxy models, and have found a best-fitting model, as well 
as a model that is best-fitting after key parameters have 
been fixed at convenient values. These models will provide a 
suitable starting point for producing fully dynamical Galaxy 
models. 

We have also shown that the main features of our mod- 
els are unchanged when we consider different disc scale- 
heights, except to the extent that they alter our prior prob- 
ability distribution on the relative contributions of the thin 
and thick discs, which alters their relative contributions in 
the models. It should therefore be noted that the (already 
weak) constraints we have on the relative contributions of 
the two discs are even weaker when the uncertainty in disc 
scale-heights is taken into account. This constraint on the ra- 
tio of thi ck to thin di sc contributions is almost entirely that 
from our lJuric et al. 

I |2008) prior, and a different choice of 
prior could result in a significantly different ratio. 

The kinematic data we consider here does not help us 
to constrain the vertical density profile of the Galactic discs, 
but it is clear that kinematic data can be used in com- 
bination with star counts to provide grea ter insight int o 
the Galactic potential above the plane (e.g. lBurnettll201(t ). 
which could then be used to improve these models. 

Applying the need for consistancy in our model allows 
us to find tighter constraints on individual parameters than 
those we begin with. This is particularly noticable in the 
constraints our posterior pdf place on the thin disc scale- 
length (-Rd,thin), the Solar radius (-Ro), and the circular ve- 
locity at the Solar radius (wo)- We find that the Galaxy's 
dark-matter halo concentration, c„/, is larger than the aver- 
age value predicted by simulations. The Galaxy's halo is less 
massive than the expected value from cosmological simula- 
tions, given the Galaxy's stellar mass (or, equivalently, the 
stellar mass is higher than would be expected). In contrast, 
the stellar mass of the Milky Way is lower than the bary- 
onic TuUy-Fisher relation would suggest given its circular 
velocity - the discrepency between the two expectations for 
the stellar mass is related to the high concentration of the 
halo. In addition to the uncertainty on individual parame- 
ters, we are able to find the correlations between different 
parameters demanded by our constraints. 

The results described in this paper are all dependent on 
the choice of parametrized model and applied constraints. In 



particular we have not attempted to take account of the ef- 
fect on the CDM halo of baryonic processes, or to consider a 
CDM halo which is not spherically symmetric. The system- 
atic uncertainties on the quantities we describe are almost 
undoubtedly larger than the quoted statistical uncertainties. 
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APPENDIX A: MODELS WITH DIFFERENT 
DISC SCALE HEIGHTS 

Por completeness we now provide a table of the mean and 
standard deviation of the parameters of our models in cases 
where we consider disc scale-heights which differ from our 
default values. These results are discussed in Section [5. II 
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Table Al. Mean and standard deviation of the parameters (upper) and derived properties (lower) of our models, for various values of 
the thin and thick disc scale- heights (/ij thin & thick respectively). This is very similar to Table O The results with our standard scale- 
heights (/id, thin = 0.3 kpc, thick = 0.9 kpc) are included in this table as well as in Table|2] Again, distances are quoted in units of kpc, 
velocities in kms""*^, masses in 10^ M0, surface densities in M0 pc"'^, densities in M0 pc~^, and -^"2,1.1 in units of (27rG) X M0 pc"'^. 



